Modeling catalytic combustion om methane using RMG-Cat

Based somewhat on the notebook to accompany the manuscript:
"Automatic generation of microkinetic mechanisms for heterogeneous catalysis" by
C. Franklin Goldsmith, School of Engineering, Brown University, franklin_goldsmith@brown.edu, and
Richard H. West, Department of Chemical Engineering, Northeastern University, r.west@northeastern.edu

As modified to prepare data for the AIChE conference presentation: 304c Incorporation of Linear Scaling Relations into Automatic Mechanism Generation for Heterogeneous Catalysis Richard H. West, Department of Chemical Engineering, Northeastern University, Boston, MA and C. Franklin Goldsmith, Engineering, Brown University, Providence, RI Tuesday, October 31, 2017 https://aiche.confex.com/aiche/2017/meetingapp.cgi/Paper/500172 https://www.slideshare.net/richardhwest/incorporation-of-linear-scaling-relations-into-automatic-mechanism-generation-for-heterogeneous-catalysis

Then further updated to model catalytic combustion of methanol.

First, we print what git commit we were on when we ran this notebook, for both the source code (RMG-Py) and the database.

In [1]:
%%bash
cd $RMGpy
pwd
git log -n1 --pretty=oneline
cd ../RMG-database
pwd
git log -n1 --pretty=oneline
/Users/rwest/Code/Cat/RMG-Py
4f7568e85b4444cceca945b6d3116ed1105070e6 Separate the Linear Scaling Relations into new function. Apply always.
/Users/rwest/Code/Cat/RMG-database
d5bf54520fc72d16d88674014df27a8fd249c3c4 Change comments (but not values) in the recommended families list.

Model generation

We start with a base input file to generate a mechanism for CH3OH plus O2. First we print the input file we'll use to generate the model.

In [2]:
%cat base/input.py
# Data sources
database(
    thermoLibraries=['surfaceThermo', 'primaryThermoLibrary', 'thermo_DFT_CCSDTF12_BAC','DFT_QCI_thermo'],
    reactionLibraries = [('Deutschmann_Ni', False)],
    seedMechanisms = [],
    kineticsDepositories = ['training'],
    kineticsFamilies = 'default',
    kineticsEstimator = 'rate rules',
    bindingEnergies = { # default values for Ni(111)
                       'C':(-5.997, 'eV/molecule'),
                       'H':(-2.778, 'eV/molecule'),
                       'O':(-4.485, 'eV/molecule'),
                       }
)

# List of species
species(
    label='X',
    reactive=True,
    structure=adjacencyList("1 X u0"),
)

species(
    label='CH4',
    reactive=True,
    structure=SMILES("[CH4]"),
)
species(
   label='O2',
   reactive=True,
   structure=adjacencyList(
       """
1 O u1 p2 c0 {2,S}
2 O u1 p2 c0 {1,S}
"""),
)

species(
    label='N2',
    reactive=False,
    structure=SMILES("N#N"),
)

species(
    label='CO2',
    reactive=True,
    structure=SMILES("O=C=O"),
)

species(
    label='H2O',
    reactive=True,
    structure=SMILES("O"),
)

species(
    label='H2',
    reactive=True,
    structure=SMILES("[H][H]"),
)

species(
    label='CO',
    reactive=True,
    structure=SMILES("[C-]#[O+]"),
)

species(
    label='C2H6',
    reactive=True,
    structure=SMILES("CC"),
)

species(
    label='CH2O',
    reactive=True,
    structure=SMILES("C=O"),
)

species(
    label='CH3',
    reactive=True,
    structure=SMILES("[CH3]"),
)

species(
    label='C3H8',
    reactive=True,
    structure=SMILES("CCC"),
)

species(
    label='H',
    reactive=True,
    structure=SMILES("[H]"),
)

species(
    label='C2H5',
    reactive=True,
    structure=SMILES("C[CH2]"),
)

species(
    label='CH3OH',
    reactive=True,
    structure=SMILES("CO"),
)

species(
    label='HCO',
    reactive=True,
    structure=SMILES("[CH]=O"),
)

species(
    label='CH3CHO',
    reactive=True,
    structure=SMILES("CC=O"),
)

species(
    label='OH',
    reactive=True,
    structure=SMILES("[OH]"),
)

species(
    label='C2H4',
    reactive=True,
    structure=SMILES("C=C"),
)

#----------
# Reaction systems
surfaceReactor(
    temperature=(1200,'K'),
    initialPressure=(1.01325, 'bar'),
    initialGasMoleFractions={
        "CH3OH": 0.11764706,
        "O2": 0.17647059,
        "N2": 0.70588235,
    },
    initialSurfaceCoverages={
        "X": 1.0,
    },
    surfaceVolumeRatio=(1.e5, 'm^-1'),
    surfaceSiteDensity=(2.9e-9, 'mol/cm^2'),
#    terminationConversion = { "CH4":0.9,},
    terminationTime=(1.0, 's'),
)

simulator(
    atol=1e-18,
    rtol=1e-12,
)

model(
    toleranceKeepInEdge=0.0,
    toleranceMoveToCore=1e-2,
    toleranceInterruptSimulation=0.1,
    maximumEdgeSpecies=100000
)

options(
    units='si',
    saveRestartPeriod=None,
    generateOutputHTML=True,
    generatePlots=False,
    saveEdgeSpecies=True,
    saveSimulationProfiles=True,
)

Then we try running it. This will take about four minutes.

In [3]:
%%bash
# python $RMGpy/rmg.py base/input.py > /dev/null
tail -n12 base/RMG.log
Saving current model edge to HTML file...
Updating RMG execution statistics...
    Execution time (DD:HH:MM:SS): 00:00:02:43
    Memory used: memory usage was unable to be logged


MODEL GENERATION COMPLETED

The final model core has 70 species and 218 reactions
The final model edge has 516 species and 977 reactions

RMG execution terminated at Wed Nov 22 15:02:38 2017

There are 70 species and 218 reactions (?)

Data processing

Next we will import some libraries and set things up to start importing and analyzing the simulation results.

In [4]:
%config InlineBackend.figure_format = 'retina'
In [5]:
%matplotlib inline
from matplotlib import pyplot as plt
import matplotlib

# The default output Type 3 (Type3) fonts can't be edited in Adobe Illustrator
# but Type 42 (TrueType) fonts can be, making it easier to move labels around
# slightly to improve layout before publication.
matplotlib.rcParams['pdf.fonttype'] = 42 

# Seaborn helps make matplotlib graphics nicer
import seaborn as sns
sns.set_style("ticks")
sns.set_context("paper", font_scale=1.5, rc={"lines.linewidth": 2.0})

import os
import re
import pandas as pd
import numpy as np
import shutil
import subprocess
import multiprocessing
In [6]:
def get_last_csv_file(job_directory):
    """
    Find the CSV file from the largest model in the provided job directory.
    
    For CSV files named `simulation_1_13.csv` you want 13 to be the highest number.
    """
    solver_directory = os.path.join(job_directory,'solver')
    csv_files = sorted([f for f in os.listdir(solver_directory) if f.endswith('.csv') ],
                       key=lambda f: int(f[:-4].split('_')[2]))
    return os.path.join(solver_directory, csv_files[-1])
    
job_directory = 'base'
get_last_csv_file(job_directory)
Out[6]:
'base/solver/simulation_1_70.csv'

We will use Pandas to import the csv file

In [7]:
last_csv_file = get_last_csv_file(job_directory)
data = pd.read_csv(last_csv_file)
data
Out[7]:
Time (s) Volume (m^3) N2 Ar He Ne X(1) CH4(2) O2(3) CO2(4) ... CH3O(51) SX(146) SX(276) SX(198) SX(250) SX(303) SX(123) SX(233) SX(87) SX(205)
0 2.659587e-24 0.098469 0.705882 0.0 0.0 0.0 0.285560 5.550766e-52 1.764706e-01 -3.000770e-120 ... 2.898438e-36 3.045506e-149 -8.214884e-112 -1.045643e-91 -1.256345e-104 -3.779405e-200 2.525411e-180 -1.666835e-229 -7.558573e-123 5.532969e-159
1 7.446844e-24 0.098469 0.705882 0.0 0.0 0.0 0.285560 8.388248e-51 1.764706e-01 5.436361e-118 ... 1.750657e-35 4.803049e-147 4.257108e-110 1.033162e-90 6.510619e-103 7.608577e-197 1.773880e-177 1.937888e-226 1.395964e-120 1.147374e-156
2 1.702136e-23 0.098469 0.705882 0.0 0.0 0.0 0.285560 9.023670e-50 1.764706e-01 2.294912e-115 ... 8.428658e-35 2.294646e-144 4.468008e-108 5.698351e-89 6.833160e-101 6.623433e-193 1.022042e-173 5.535459e-221 5.818797e-118 3.139309e-153
3 3.617039e-23 0.098469 0.705882 0.0 0.0 0.0 0.285560 8.294299e-49 1.764706e-01 7.381584e-113 ... 3.681017e-34 7.792782e-142 3.567230e-106 2.242674e-87 5.455553e-99 3.851828e-189 3.411043e-170 3.260231e-216 1.872945e-115 5.010657e-150
4 7.446844e-23 0.098469 0.705882 0.0 0.0 0.0 0.285560 7.095944e-48 1.764706e-01 2.111530e-110 ... 1.536752e-33 2.286727e-139 2.540984e-104 7.923873e-86 3.886061e-97 1.869238e-185 8.816793e-167 1.417965e-211 5.359923e-113 6.326241e-147
5 1.510646e-22 0.098469 0.705882 0.0 0.0 0.0 0.285560 5.867062e-47 1.764706e-01 5.709668e-108 ... 6.278133e-33 6.260534e-137 1.714242e-102 2.662090e-84 2.621680e-95 8.322394e-182 2.022951e-163 5.338346e-207 1.449660e-110 7.172265e-144
6 3.042568e-22 0.098469 0.705882 0.0 0.0 0.0 0.285560 4.770997e-46 1.764706e-01 1.501952e-105 ... 2.537722e-32 1.656961e-134 1.126187e-100 8.726700e-83 1.722337e-93 3.552828e-178 4.382233e-160 1.873706e-202 3.813806e-108 7.722493e-141
7 6.106412e-22 0.098469 0.705882 0.0 0.0 0.0 0.285560 3.847982e-45 1.764706e-01 3.897430e-103 ... 1.020407e-31 4.312744e-132 7.302108e-99 2.826288e-81 1.116750e-91 1.485529e-174 9.228764e-157 6.353288e-198 9.897021e-106 8.107436e-138
8 1.223410e-21 0.098469 0.705882 0.0 0.0 0.0 0.285560 3.090908e-44 1.764706e-01 1.004510e-100 ... 4.092288e-31 1.113236e-129 4.703832e-97 9.098482e-80 7.193817e-90 6.147598e-171 1.916526e-153 2.117645e-193 2.550894e-103 8.405781e-135
9 2.448948e-21 0.098469 0.705882 0.0 0.0 0.0 0.285560 2.477745e-43 1.764706e-01 2.580245e-98 ... 1.639050e-30 2.861687e-127 3.020246e-95 2.920242e-78 4.619020e-88 2.531019e-167 3.952399e-150 6.998433e-189 6.552469e-101 8.661052e-132
10 4.900023e-21 0.098469 0.705882 0.0 0.0 0.0 0.285560 1.984206e-42 1.764706e-01 6.616585e-96 ... 6.560470e-30 7.341065e-125 1.936097e-93 9.358763e-77 2.960974e-86 1.039369e-163 8.122648e-147 2.303023e-184 1.680277e-98 8.896429e-129
11 9.802174e-21 0.098469 0.705882 0.0 0.0 0.0 0.285560 1.588169e-41 1.764706e-01 1.695275e-93 ... 2.625042e-29 1.881254e-122 1.240108e-91 2.997044e-75 1.896560e-84 4.262716e-160 1.666405e-143 7.562605e-180 4.305152e-96 9.124053e-126
12 1.960648e-20 0.098469 0.705882 0.0 0.0 0.0 0.285560 1.270857e-40 1.764706e-01 4.341730e-91 ... 1.050188e-28 4.818497e-120 7.939911e-90 9.594127e-74 1.214290e-82 1.747127e-156 3.415755e-140 2.480746e-175 1.102585e-93 9.350254e-123
13 3.921508e-20 0.098469 0.705882 0.0 0.0 0.0 0.285560 1.016814e-39 1.764706e-01 1.111715e-88 ... 4.201093e-28 1.233854e-117 5.082574e-88 3.070694e-72 7.773015e-81 7.158511e-153 6.998485e-137 8.133205e-171 2.823209e-91 9.578346e-120
14 7.843229e-20 0.098469 0.705882 0.0 0.0 0.0 0.285560 8.135029e-39 1.764706e-01 2.846278e-86 ... 1.680506e-27 3.159073e-115 3.253177e-86 9.827139e-71 4.975213e-79 2.932584e-149 1.433595e-133 2.665780e-166 7.228150e-89 9.810084e-117
15 1.568667e-19 0.098469 0.705882 0.0 0.0 0.0 0.285560 6.508229e-38 1.764706e-01 7.286794e-84 ... 6.722160e-27 8.087747e-113 2.082139e-84 3.144831e-69 3.184271e-77 1.201272e-145 2.936298e-130 8.736280e-162 1.850489e-86 1.004642e-113
16 3.137355e-19 0.098469 0.705882 0.0 0.0 0.0 0.285560 5.206665e-37 1.764706e-01 1.865436e-81 ... 2.688891e-26 2.070530e-110 1.332603e-82 1.006369e-67 2.037951e-75 4.920526e-142 6.013773e-127 2.862822e-157 4.737298e-84 1.028786e-110
17 6.274732e-19 0.098469 0.705882 0.0 0.0 0.0 0.285560 4.165365e-36 1.764706e-01 4.775417e-79 ... 1.075562e-25 5.300640e-108 8.528768e-81 3.220418e-66 1.304262e-73 2.015420e-138 1.231617e-123 9.380737e-153 1.212724e-81 1.053469e-107
18 1.254949e-18 0.098469 0.705882 0.0 0.0 0.0 0.285560 3.332305e-35 1.764706e-01 1.222431e-76 ... 4.302259e-25 1.356974e-105 5.458448e-79 1.030539e-64 8.346771e-72 8.254685e-135 2.522232e-120 3.073622e-148 3.104387e-79 1.078696e-104
19 2.509899e-18 0.098469 0.705882 0.0 0.0 0.0 0.285560 2.665848e-34 1.764706e-01 3.129004e-74 ... 1.720906e-24 3.473863e-103 3.493420e-77 3.297727e-63 5.341236e-70 3.380679e-131 2.371121e-116 1.749979e-142 7.946197e-77 2.949734e-101
20 5.019801e-18 0.098469 0.705882 0.0 0.0 0.0 0.285560 2.132680e-33 1.764706e-01 8.008066e-72 ... 6.883627e-24 8.893084e-101 2.235795e-75 1.055270e-61 3.417483e-68 1.384356e-127 4.856631e-113 5.731453e-138 2.033688e-74 3.021574e-98
21 1.003960e-17 0.098469 0.705882 0.0 0.0 0.0 0.285560 1.706143e-32 1.764706e-01 2.048943e-69 ... 2.753452e-23 2.276620e-98 1.430914e-73 3.376843e-60 2.186023e-66 5.667268e-124 9.941107e-110 1.876083e-133 5.203471e-72 3.092446e-95
22 2.007921e-17 0.098469 0.705882 0.0 0.0 0.0 0.285560 1.364912e-31 1.764706e-01 5.239548e-67 ... 1.101381e-22 5.828091e-96 9.157898e-72 1.080574e-58 1.397562e-64 2.318811e-120 2.033778e-106 6.134457e-129 1.330671e-69 3.163290e-92
23 3.011881e-17 0.098469 0.705882 0.0 0.0 0.0 0.285560 4.052074e-31 1.764706e-01 3.977051e-66 ... 2.340434e-22 4.471629e-95 5.226004e-71 5.169632e-58 7.972531e-64 2.679590e-119 1.134793e-105 3.227929e-128 1.007226e-68 1.901276e-91
24 4.015842e-17 0.098469 0.705882 0.0 0.0 0.0 0.285560 8.829240e-31 1.764706e-01 1.956381e-65 ... 3.992506e-22 2.347741e-94 2.044560e-70 1.710962e-57 3.117883e-63 2.119558e-118 1.528775e-104 9.249952e-127 5.070733e-68 1.896139e-90
25 6.023763e-17 0.098469 0.705882 0.0 0.0 0.0 0.285560 2.562846e-30 1.764706e-01 2.482671e-64 ... 8.265977e-22 3.984892e-93 1.712316e-69 9.982850e-57 2.608369e-62 9.799278e-117 1.263060e-102 3.395889e-124 6.953826e-67 9.182620e-89
26 8.031684e-17 0.098469 0.705882 0.0 0.0 0.0 0.285560 5.604638e-30 1.764706e-01 1.269641e-63 ... 1.405826e-21 2.142298e-92 6.729606e-69 3.321393e-56 1.024531e-61 8.395371e-116 1.008908e-101 4.479775e-123 3.416643e-66 6.399979e-88
27 1.003961e-16 0.098469 0.705882 0.0 0.0 0.0 0.285560 1.051504e-29 1.764706e-01 4.975969e-63 ... 2.147785e-21 8.983934e-92 2.113200e-68 8.955630e-56 3.215223e-61 5.528779e-115 7.139219e-101 5.394558e-122 1.359703e-65 3.829303e-87
28 1.405545e-16 0.098469 0.705882 0.0 0.0 0.0 0.285560 2.877088e-29 1.764706e-01 1.165412e-61 ... 4.130114e-21 2.394535e-90 1.923599e-67 5.279505e-55 2.918093e-60 1.957278e-112 6.496645e-99 4.698128e-119 3.163474e-64 2.193521e-85
29 1.807129e-16 0.098469 0.705882 0.0 0.0 0.0 0.285560 6.065219e-29 1.764706e-01 7.022635e-61 ... 6.773376e-21 1.448145e-89 7.984640e-67 1.780673e-54 1.209276e-59 2.013615e-111 5.825554e-98 5.484548e-118 1.877377e-63 1.771828e-84
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
11425 9.600000e-01 0.098469 0.705882 0.0 0.0 0.0 0.240260 2.124707e-10 5.082856e-13 9.997758e-02 ... 4.290758e-12 4.742187e-15 6.848915e-12 2.622838e-14 1.384592e-10 3.404114e-17 1.780583e-09 5.873223e-13 1.694978e-09 6.718834e-11
11426 9.614785e-01 0.098469 0.705882 0.0 0.0 0.0 0.240280 2.118370e-10 5.078101e-13 9.998590e-02 ... 4.284520e-12 4.735093e-15 6.806761e-12 2.621597e-14 1.375960e-10 3.380244e-17 1.780510e-09 5.840538e-13 1.696645e-09 6.722294e-11
11427 9.629569e-01 0.098469 0.705882 0.0 0.0 0.0 0.240300 2.112072e-10 5.073359e-13 9.999420e-02 ... 4.278311e-12 4.728016e-15 6.764868e-12 2.620358e-14 1.367381e-10 3.356545e-17 1.780436e-09 5.808029e-13 1.698311e-09 6.725752e-11
11428 9.644354e-01 0.098469 0.705882 0.0 0.0 0.0 0.240320 2.105811e-10 5.068628e-13 1.000025e-01 ... 4.272132e-12 4.720959e-15 6.723235e-12 2.619121e-14 1.358855e-10 3.333016e-17 1.780361e-09 5.775696e-13 1.699976e-09 6.729206e-11
11429 9.659139e-01 0.098469 0.705882 0.0 0.0 0.0 0.240340 2.099589e-10 5.063909e-13 1.000108e-01 ... 4.265982e-12 4.713919e-15 6.681860e-12 2.617888e-14 1.350382e-10 3.309657e-17 1.780285e-09 5.743536e-13 1.701640e-09 6.732656e-11
11430 9.673923e-01 0.098469 0.705882 0.0 0.0 0.0 0.240360 2.093404e-10 5.059201e-13 1.000190e-01 ... 4.259862e-12 4.706898e-15 6.640742e-12 2.616656e-14 1.341960e-10 3.286465e-17 1.780208e-09 5.711551e-13 1.703303e-09 6.736103e-11
11431 9.688708e-01 0.098469 0.705882 0.0 0.0 0.0 0.240380 2.087256e-10 5.054506e-13 1.000273e-01 ... 4.253770e-12 4.699896e-15 6.599879e-12 2.615427e-14 1.333591e-10 3.263440e-17 1.780130e-09 5.679737e-13 1.704965e-09 6.739546e-11
11432 9.703492e-01 0.098469 0.705882 0.0 0.0 0.0 0.240399 2.081145e-10 5.049822e-13 1.000355e-01 ... 4.247708e-12 4.692911e-15 6.559269e-12 2.614201e-14 1.325273e-10 3.240580e-17 1.780051e-09 5.648096e-13 1.706626e-09 6.742986e-11
11433 9.718277e-01 0.098469 0.705882 0.0 0.0 0.0 0.240419 2.075071e-10 5.045150e-13 1.000437e-01 ... 4.241675e-12 4.685944e-15 6.518911e-12 2.612977e-14 1.317007e-10 3.217885e-17 1.779971e-09 5.616625e-13 1.708286e-09 6.746423e-11
11434 9.733062e-01 0.098469 0.705882 0.0 0.0 0.0 0.240439 2.069034e-10 5.040489e-13 1.000519e-01 ... 4.235670e-12 4.678996e-15 6.478803e-12 2.611755e-14 1.308791e-10 3.195352e-17 1.779890e-09 5.585324e-13 1.709946e-09 6.749856e-11
11435 9.747846e-01 0.098469 0.705882 0.0 0.0 0.0 0.240459 2.063034e-10 5.035840e-13 1.000601e-01 ... 4.229694e-12 4.672065e-15 6.438944e-12 2.610536e-14 1.300626e-10 3.172980e-17 1.779808e-09 5.554192e-13 1.711604e-09 6.753286e-11
11436 9.762631e-01 0.098469 0.705882 0.0 0.0 0.0 0.240478 2.057069e-10 5.031203e-13 1.000683e-01 ... 4.223746e-12 4.665152e-15 6.399332e-12 2.609319e-14 1.292511e-10 3.150770e-17 1.779726e-09 5.523229e-13 1.713261e-09 6.756713e-11
11437 9.777415e-01 0.098469 0.705882 0.0 0.0 0.0 0.240498 2.051141e-10 5.026577e-13 1.000764e-01 ... 4.217827e-12 4.658257e-15 6.359966e-12 2.608104e-14 1.284447e-10 3.128718e-17 1.779642e-09 5.492432e-13 1.714918e-09 6.760136e-11
11438 9.792200e-01 0.098469 0.705882 0.0 0.0 0.0 0.240517 2.045249e-10 5.021962e-13 1.000846e-01 ... 4.211935e-12 4.651380e-15 6.320843e-12 2.606892e-14 1.276432e-10 3.106825e-17 1.779557e-09 5.461803e-13 1.716574e-09 6.763556e-11
11439 9.806985e-01 0.098469 0.705882 0.0 0.0 0.0 0.240537 2.039392e-10 5.017358e-13 1.000927e-01 ... 4.206072e-12 4.644520e-15 6.281963e-12 2.605683e-14 1.268466e-10 3.085088e-17 1.779471e-09 5.431338e-13 1.718228e-09 6.766972e-11
11440 9.821769e-01 0.098469 0.705882 0.0 0.0 0.0 0.240557 2.033571e-10 5.012766e-13 1.001008e-01 ... 4.200236e-12 4.637678e-15 6.243324e-12 2.604475e-14 1.260550e-10 3.063507e-17 1.779385e-09 5.401039e-13 1.719882e-09 6.770385e-11
11441 9.836554e-01 0.098469 0.705882 0.0 0.0 0.0 0.240576 2.027784e-10 5.008186e-13 1.001090e-01 ... 4.194428e-12 4.630853e-15 6.204924e-12 2.603270e-14 1.252682e-10 3.042081e-17 1.779297e-09 5.370904e-13 1.721535e-09 6.773795e-11
11442 9.851338e-01 0.098469 0.705882 0.0 0.0 0.0 0.240595 2.022034e-10 5.003616e-13 1.001170e-01 ... 4.188648e-12 4.624046e-15 6.166762e-12 2.602067e-14 1.244863e-10 3.020808e-17 1.779208e-09 5.340931e-13 1.723187e-09 6.777202e-11
11443 9.866123e-01 0.098469 0.705882 0.0 0.0 0.0 0.240615 2.016317e-10 4.999058e-13 1.001251e-01 ... 4.182895e-12 4.617256e-15 6.128837e-12 2.600867e-14 1.237092e-10 2.999687e-17 1.779119e-09 5.311121e-13 1.724839e-09 6.780605e-11
11444 9.880908e-01 0.098469 0.705882 0.0 0.0 0.0 0.240634 2.010636e-10 4.994510e-13 1.001332e-01 ... 4.177169e-12 4.610483e-15 6.091146e-12 2.599669e-14 1.229369e-10 2.978717e-17 1.779028e-09 5.281472e-13 1.726489e-09 6.784005e-11
11445 9.895692e-01 0.098469 0.705882 0.0 0.0 0.0 0.240654 2.004989e-10 4.989974e-13 1.001412e-01 ... 4.171470e-12 4.603727e-15 6.053689e-12 2.598473e-14 1.221693e-10 2.957898e-17 1.778937e-09 5.251984e-13 1.728138e-09 6.787402e-11
11446 9.910477e-01 0.098469 0.705882 0.0 0.0 0.0 0.240673 1.999376e-10 4.985448e-13 1.001493e-01 ... 4.165798e-12 4.596989e-15 6.016464e-12 2.597280e-14 1.214065e-10 2.937227e-17 1.778845e-09 5.222656e-13 1.729787e-09 6.790795e-11
11447 9.925261e-01 0.098469 0.705882 0.0 0.0 0.0 0.240692 1.993797e-10 4.980934e-13 1.001573e-01 ... 4.160153e-12 4.590267e-15 5.979469e-12 2.596089e-14 1.206484e-10 2.916704e-17 1.778752e-09 5.193486e-13 1.731434e-09 6.794185e-11
11448 9.938568e-01 0.098469 0.705882 0.0 0.0 0.0 0.240710 1.988806e-10 4.976880e-13 1.001645e-01 ... 4.155095e-12 4.584233e-15 5.946370e-12 2.595018e-14 1.199701e-10 2.898359e-17 1.778667e-09 5.167369e-13 1.732917e-09 6.797234e-11
11449 9.951874e-01 0.098469 0.705882 0.0 0.0 0.0 0.240727 1.983841e-10 4.972836e-13 1.001717e-01 ... 4.150058e-12 4.578211e-15 5.913456e-12 2.593950e-14 1.192956e-10 2.880132e-17 1.778581e-09 5.141379e-13 1.734398e-09 6.800279e-11
11450 9.965180e-01 0.098469 0.705882 0.0 0.0 0.0 0.240744 1.978904e-10 4.968799e-13 1.001789e-01 ... 4.145043e-12 4.572204e-15 5.880725e-12 2.592884e-14 1.186248e-10 2.862023e-17 1.778495e-09 5.115516e-13 1.735879e-09 6.803323e-11
11451 9.978486e-01 0.098469 0.705882 0.0 0.0 0.0 0.240762 1.973993e-10 4.964772e-13 1.001861e-01 ... 4.140050e-12 4.566210e-15 5.848176e-12 2.591819e-14 1.179578e-10 2.844030e-17 1.778408e-09 5.089780e-13 1.737359e-09 6.806363e-11
11452 9.991792e-01 0.098469 0.705882 0.0 0.0 0.0 0.240779 1.969110e-10 4.960753e-13 1.001933e-01 ... 4.135077e-12 4.560230e-15 5.815809e-12 2.590756e-14 1.172945e-10 2.826152e-17 1.778321e-09 5.064169e-13 1.738838e-09 6.809401e-11
11453 1.000000e+00 0.098469 0.705882 0.0 0.0 0.0 0.240790 1.966111e-10 4.958279e-13 1.001977e-01 ... 4.132020e-12 4.556547e-15 5.795933e-12 2.590102e-14 1.168871e-10 2.815182e-17 1.778266e-09 5.048434e-13 1.739750e-09 6.811274e-11
11454 1.000377e+00 0.098469 0.705882 0.0 0.0 0.0 0.240794 1.964737e-10 4.957144e-13 1.001997e-01 ... 4.130619e-12 4.554859e-15 5.786833e-12 2.589801e-14 1.167006e-10 2.810161e-17 1.778241e-09 5.041226e-13 1.740169e-09 6.812133e-11

11455 rows × 72 columns

In [8]:
def get_pandas_data(job_directory):
    """
    Get the last CSV file from the provided job directory,
    import it into a Pandas data frame, and tidy it up a bit.
    """
    last_csv_file = get_last_csv_file(job_directory)
    data = pd.read_csv(last_csv_file)
    
    # Make the Time into an index and remove it as a column
    data.index = data['Time (s)']
    del data['Time (s)']
    # Remove inerts that RMG added automatically but we're not using
    for i in 'Ar He Ne'.split():
        del data[i]
    # Remove the Volume column
    del data['Volume (m^3)']
    # Set any amounts that are slightly negative (within the ODE solver's ATOL tolerance), equal to zero
    # to allow 'area' plots without warnings or errors.
    # Anything more negative than -1e-16 probably indicates a bug in RMG and should not be hidden here ;-)
    data[(data<0) & (data>-1e-16)] = 0
    return data
In [9]:
def rename_columns(data):
    """
    Removes the number (##) from the end of the column names, in place,
    unless doing so would make the names collide.
    Also renames a few species so the plot labels match the names in the manuscript.
    """
    import re
    old = data.columns
    new = [re.sub('\(\d+\)$','',n) for n in old]
    # don't translate names that would no longer be unique
    mapping = {k:v for k,v in zip(old,new) if new.count(v)==1}
    data.rename(columns=mapping, inplace=True)
    
    # Now change a few species that are named differently in the manuscript
    # than in the thermodynamics database used to build the model,
    # so that the plot labels match the manuscript.
    mapping = {'COX':'COvdwX', 'OCX': 'CO=X', 'C2H3X':'CH3CX', 'C2H3OX':'CH3CXO'}
    data.rename(columns=mapping, inplace=True)
In [10]:
data1 = get_pandas_data('base')
rename_columns(data1)
data1.columns
Out[10]:
Index([u'N2', u'X', u'CH4', u'O2', u'CO2', u'H2O', u'H2', u'CO', u'C2H6',
       u'CH2O', u'CH3', u'C3H8', u'H', u'C2H5', u'CH3OH', u'HCO', u'CH3CHO',
       u'OH', u'C2H4', u'CH4OX', u'OX', u'HX', u'CH3OX(39)', u'HOX', u'OCCO',
       u'CH2OX', u'CHOX', u'H2X', u'H2OX', u'CH3X', u'CH2X', u'CHX', u'CXHO',
       u'CO=X', u'COvdwX', u'CH3OX(40)', u'HO2X', u'CH4X', u'CH2O2X',
       u'CH3O2X(72)', u'CH3O2X(89)', u'CX', u'SX(68)', u'CO2X', u'SX(114)',
       u'SX(93)', u'SX(141)', u'C2H4OX', u'SX(166)', u'SX(138)', u'CH3O2X(61)',
       u'SX(169)', u'SX(98)', u'C2H5OX', u'SX(92)', u'SX(224)', u'HOCXO',
       u'CH3O', u'SX(146)', u'SX(276)', u'SX(198)', u'SX(250)', u'SX(303)',
       u'SX(123)', u'SX(233)', u'SX(87)', u'SX(205)'],
      dtype='object')
In [11]:
# Test it with some plots
data1[['CH3OH', 'O2']].plot.line()
data1[['CO', 'H2']].plot.line()
data1[['H2O']].plot.line()
Out[11]:
<matplotlib.axes._subplots.AxesSubplot at 0x11b85b6d0>
In [12]:
species_names = data1.columns
# just the gas phase species that aren't always zero:
gas_phase = [n for n in species_names if 'X' not in n and (data1[n]>0).any()]
# all the surface species
surface_phase = [n for n in species_names if 'X' in n]
surface_phase.remove('X')
data1[gas_phase].plot.line()
data1[surface_phase].plot.line()
Out[12]:
<matplotlib.axes._subplots.AxesSubplot at 0x11c749990>
In [13]:
print "Significant species (those that exceed 0.001 mol at some point)"
significant = [n for n in data1.columns if(data1[n]>0.001).any()]
with sns.color_palette("hls", len(significant)):
    data1[significant].plot.area(legend='reverse')
Significant species (those that exceed 0.001 mol at some point)
In [14]:
lim = 1e-4
surface = [n for n in data1.columns if 'X' in n and n!='X' and (data1[n]>lim).any() ]
print "The {} surface species that exceed {:.1e} mol at some point".format(len(surface), lim)
total_sites = max(data1['X'])
with sns.color_palette('Set3',len(surface)):
    (data1[surface]/total_sites).plot.area(legend='reverse')
    plt.ylabel('site fraction')
    plt.xlim(0,1e-2) ## notice this is much less than 1 seconde
The 12 surface species that exceed 1.0e-04 mol at some point

Effect of binding energies

Now we will use that template 'base' input file to create a ton of other input files with different binding energies, then run them all in RMG-Cat, then process the results.

In [15]:
# First, make a series of input files in separate directories

with open(os.path.join('base', 'input.py')) as infile:
    input_file = infile.read()
    
base_directory = 'binding_energies'
def directory(carbon,oxygen):
    return os.path.join(base_directory, "c{:.3f}o{:.3f}".format(carbon,oxygen))

def make_input(binding_energies):
    """
    Make an input file for the given (carbon,oxygen) tuple (or iterable) of binding energies
    and return the name of the directory in which it is saved.
    """
    carbon, oxygen = binding_energies
    output = input_file
    out_dir = directory(carbon, oxygen)
    carbon_string = "'C':({:f}, 'eV/molecule')".format(carbon)
    output = re.sub("'C':\(.*?, 'eV/molecule'\)", carbon_string, output)
    oxygen_string = "'O':({:f}, 'eV/molecule')".format(oxygen)
    output = re.sub("'O':\(.*?, 'eV/molecule'\)", oxygen_string, output)
    os.path.exists(out_dir) or os.makedirs(out_dir)
    out_file = os.path.join(out_dir, 'input.py')
    with open(out_file,'w') as outfile:
        outfile.write(output)
    shutil.copy(os.path.join('base','run.sh'), out_dir)
    return out_dir

    
print make_input((-8,-3.5))
    
binding_energies/c-8.000o-3.500
In [16]:
def run_simulation(binding_energies):
    """
    Assuming a job file already exists, run it.  This one is local.
    Takes a tuple of binding energies, (carbon, oxygen)
    """
    carbon, oxygen = binding_energies
    job_directory = directory(carbon, oxygen)
    print job_directory
    assert os.path.exists(job_directory)
    return subprocess.check_call('./run.sh', cwd=job_directory)
# a test experiment = (-8,-3.5) make_input(experiment) run_simulation(experiment)
In [17]:
# Revised range
carbon_range = (-7.5, -2)
oxygen_range = (-6.5, -1.5)
grid_size = 9
mesh  = np.mgrid[carbon_range[0]:carbon_range[1]:grid_size*1j, 
                 oxygen_range[0]:oxygen_range[1]:grid_size*1j]

with sns.axes_style("whitegrid"):
    plt.axis('square')
    plt.xlim(carbon_range)
    plt.ylim(oxygen_range)
plt.show()
    
mesh
Out[17]:
array([[[-7.5   , -7.5   , -7.5   , -7.5   , -7.5   , -7.5   , -7.5   ,
         -7.5   , -7.5   ],
        [-6.8125, -6.8125, -6.8125, -6.8125, -6.8125, -6.8125, -6.8125,
         -6.8125, -6.8125],
        [-6.125 , -6.125 , -6.125 , -6.125 , -6.125 , -6.125 , -6.125 ,
         -6.125 , -6.125 ],
        [-5.4375, -5.4375, -5.4375, -5.4375, -5.4375, -5.4375, -5.4375,
         -5.4375, -5.4375],
        [-4.75  , -4.75  , -4.75  , -4.75  , -4.75  , -4.75  , -4.75  ,
         -4.75  , -4.75  ],
        [-4.0625, -4.0625, -4.0625, -4.0625, -4.0625, -4.0625, -4.0625,
         -4.0625, -4.0625],
        [-3.375 , -3.375 , -3.375 , -3.375 , -3.375 , -3.375 , -3.375 ,
         -3.375 , -3.375 ],
        [-2.6875, -2.6875, -2.6875, -2.6875, -2.6875, -2.6875, -2.6875,
         -2.6875, -2.6875],
        [-2.    , -2.    , -2.    , -2.    , -2.    , -2.    , -2.    ,
         -2.    , -2.    ]],

       [[-6.5   , -5.875 , -5.25  , -4.625 , -4.    , -3.375 , -2.75  ,
         -2.125 , -1.5   ],
        [-6.5   , -5.875 , -5.25  , -4.625 , -4.    , -3.375 , -2.75  ,
         -2.125 , -1.5   ],
        [-6.5   , -5.875 , -5.25  , -4.625 , -4.    , -3.375 , -2.75  ,
         -2.125 , -1.5   ],
        [-6.5   , -5.875 , -5.25  , -4.625 , -4.    , -3.375 , -2.75  ,
         -2.125 , -1.5   ],
        [-6.5   , -5.875 , -5.25  , -4.625 , -4.    , -3.375 , -2.75  ,
         -2.125 , -1.5   ],
        [-6.5   , -5.875 , -5.25  , -4.625 , -4.    , -3.375 , -2.75  ,
         -2.125 , -1.5   ],
        [-6.5   , -5.875 , -5.25  , -4.625 , -4.    , -3.375 , -2.75  ,
         -2.125 , -1.5   ],
        [-6.5   , -5.875 , -5.25  , -4.625 , -4.    , -3.375 , -2.75  ,
         -2.125 , -1.5   ],
        [-6.5   , -5.875 , -5.25  , -4.625 , -4.    , -3.375 , -2.75  ,
         -2.125 , -1.5   ]]])
In [18]:
experiments = mesh.reshape((2,-1)).T
experiments
Out[18]:
array([[-7.5   , -6.5   ],
       [-7.5   , -5.875 ],
       [-7.5   , -5.25  ],
       [-7.5   , -4.625 ],
       [-7.5   , -4.    ],
       [-7.5   , -3.375 ],
       [-7.5   , -2.75  ],
       [-7.5   , -2.125 ],
       [-7.5   , -1.5   ],
       [-6.8125, -6.5   ],
       [-6.8125, -5.875 ],
       [-6.8125, -5.25  ],
       [-6.8125, -4.625 ],
       [-6.8125, -4.    ],
       [-6.8125, -3.375 ],
       [-6.8125, -2.75  ],
       [-6.8125, -2.125 ],
       [-6.8125, -1.5   ],
       [-6.125 , -6.5   ],
       [-6.125 , -5.875 ],
       [-6.125 , -5.25  ],
       [-6.125 , -4.625 ],
       [-6.125 , -4.    ],
       [-6.125 , -3.375 ],
       [-6.125 , -2.75  ],
       [-6.125 , -2.125 ],
       [-6.125 , -1.5   ],
       [-5.4375, -6.5   ],
       [-5.4375, -5.875 ],
       [-5.4375, -5.25  ],
       [-5.4375, -4.625 ],
       [-5.4375, -4.    ],
       [-5.4375, -3.375 ],
       [-5.4375, -2.75  ],
       [-5.4375, -2.125 ],
       [-5.4375, -1.5   ],
       [-4.75  , -6.5   ],
       [-4.75  , -5.875 ],
       [-4.75  , -5.25  ],
       [-4.75  , -4.625 ],
       [-4.75  , -4.    ],
       [-4.75  , -3.375 ],
       [-4.75  , -2.75  ],
       [-4.75  , -2.125 ],
       [-4.75  , -1.5   ],
       [-4.0625, -6.5   ],
       [-4.0625, -5.875 ],
       [-4.0625, -5.25  ],
       [-4.0625, -4.625 ],
       [-4.0625, -4.    ],
       [-4.0625, -3.375 ],
       [-4.0625, -2.75  ],
       [-4.0625, -2.125 ],
       [-4.0625, -1.5   ],
       [-3.375 , -6.5   ],
       [-3.375 , -5.875 ],
       [-3.375 , -5.25  ],
       [-3.375 , -4.625 ],
       [-3.375 , -4.    ],
       [-3.375 , -3.375 ],
       [-3.375 , -2.75  ],
       [-3.375 , -2.125 ],
       [-3.375 , -1.5   ],
       [-2.6875, -6.5   ],
       [-2.6875, -5.875 ],
       [-2.6875, -5.25  ],
       [-2.6875, -4.625 ],
       [-2.6875, -4.    ],
       [-2.6875, -3.375 ],
       [-2.6875, -2.75  ],
       [-2.6875, -2.125 ],
       [-2.6875, -1.5   ],
       [-2.    , -6.5   ],
       [-2.    , -5.875 ],
       [-2.    , -5.25  ],
       [-2.    , -4.625 ],
       [-2.    , -4.    ],
       [-2.    , -3.375 ],
       [-2.    , -2.75  ],
       [-2.    , -2.125 ],
       [-2.    , -1.5   ]])
In [19]:
with sns.axes_style("whitegrid"):
    plt.axis('square')
    plt.xlim(carbon_range)
    plt.ylim(oxygen_range)
    plt.plot(*experiments.T, marker='o', linestyle='none')
In [20]:
map(make_input, experiments)
Out[20]:
['binding_energies/c-7.500o-6.500',
 'binding_energies/c-7.500o-5.875',
 'binding_energies/c-7.500o-5.250',
 'binding_energies/c-7.500o-4.625',
 'binding_energies/c-7.500o-4.000',
 'binding_energies/c-7.500o-3.375',
 'binding_energies/c-7.500o-2.750',
 'binding_energies/c-7.500o-2.125',
 'binding_energies/c-7.500o-1.500',
 'binding_energies/c-6.812o-6.500',
 'binding_energies/c-6.812o-5.875',
 'binding_energies/c-6.812o-5.250',
 'binding_energies/c-6.812o-4.625',
 'binding_energies/c-6.812o-4.000',
 'binding_energies/c-6.812o-3.375',
 'binding_energies/c-6.812o-2.750',
 'binding_energies/c-6.812o-2.125',
 'binding_energies/c-6.812o-1.500',
 'binding_energies/c-6.125o-6.500',
 'binding_energies/c-6.125o-5.875',
 'binding_energies/c-6.125o-5.250',
 'binding_energies/c-6.125o-4.625',
 'binding_energies/c-6.125o-4.000',
 'binding_energies/c-6.125o-3.375',
 'binding_energies/c-6.125o-2.750',
 'binding_energies/c-6.125o-2.125',
 'binding_energies/c-6.125o-1.500',
 'binding_energies/c-5.438o-6.500',
 'binding_energies/c-5.438o-5.875',
 'binding_energies/c-5.438o-5.250',
 'binding_energies/c-5.438o-4.625',
 'binding_energies/c-5.438o-4.000',
 'binding_energies/c-5.438o-3.375',
 'binding_energies/c-5.438o-2.750',
 'binding_energies/c-5.438o-2.125',
 'binding_energies/c-5.438o-1.500',
 'binding_energies/c-4.750o-6.500',
 'binding_energies/c-4.750o-5.875',
 'binding_energies/c-4.750o-5.250',
 'binding_energies/c-4.750o-4.625',
 'binding_energies/c-4.750o-4.000',
 'binding_energies/c-4.750o-3.375',
 'binding_energies/c-4.750o-2.750',
 'binding_energies/c-4.750o-2.125',
 'binding_energies/c-4.750o-1.500',
 'binding_energies/c-4.062o-6.500',
 'binding_energies/c-4.062o-5.875',
 'binding_energies/c-4.062o-5.250',
 'binding_energies/c-4.062o-4.625',
 'binding_energies/c-4.062o-4.000',
 'binding_energies/c-4.062o-3.375',
 'binding_energies/c-4.062o-2.750',
 'binding_energies/c-4.062o-2.125',
 'binding_energies/c-4.062o-1.500',
 'binding_energies/c-3.375o-6.500',
 'binding_energies/c-3.375o-5.875',
 'binding_energies/c-3.375o-5.250',
 'binding_energies/c-3.375o-4.625',
 'binding_energies/c-3.375o-4.000',
 'binding_energies/c-3.375o-3.375',
 'binding_energies/c-3.375o-2.750',
 'binding_energies/c-3.375o-2.125',
 'binding_energies/c-3.375o-1.500',
 'binding_energies/c-2.688o-6.500',
 'binding_energies/c-2.688o-5.875',
 'binding_energies/c-2.688o-5.250',
 'binding_energies/c-2.688o-4.625',
 'binding_energies/c-2.688o-4.000',
 'binding_energies/c-2.688o-3.375',
 'binding_energies/c-2.688o-2.750',
 'binding_energies/c-2.688o-2.125',
 'binding_energies/c-2.688o-1.500',
 'binding_energies/c-2.000o-6.500',
 'binding_energies/c-2.000o-5.875',
 'binding_energies/c-2.000o-5.250',
 'binding_energies/c-2.000o-4.625',
 'binding_energies/c-2.000o-4.000',
 'binding_energies/c-2.000o-3.375',
 'binding_energies/c-2.000o-2.750',
 'binding_energies/c-2.000o-2.125',
 'binding_energies/c-2.000o-1.500']
In [21]:
## Now run the simulations using a pool.
## Don't run this cell unless you have a while to wait!! ###

# pool = multiprocessing.Pool()
# result = pool.map(run_simulation, experiments)

## Instead, upload the input files to a cluster and run `start_all.sh`
## The script `stage_all.sh` will stage (add) the results in git
## so you can commit and get them back to analyze with the subsequent cells:
In [22]:
base_directory = 'binding_energies'
# base_directory = 'binding_energies_local'
In [23]:
def get_data(experiment):
    carbon, oxygen = experiment
    directory(carbon,oxygen)
    try:
        data = get_pandas_data(directory(carbon,oxygen))
    except OSError:
        return None
    rename_columns(data)
    return data
In [24]:
datas = {tuple(e): get_data(e) for e in experiments}
In [25]:
def get_max_co2(experiment):
    data = datas[tuple(experiment)]
    if data is None:
        return np.nan
    return data[['CO2']].max()
highest_co2 = max([float(get_max_co2(e)) for e in experiments])

Fit an exponential

For this first attempt to extract "rate" we fit an exponential growth curve to the normalized CO2 concentration profile of each simulation. First we plot all the curves on one plot, then we'll plot each with its fitted exponential.

In [26]:
import seaborn as sns
plt.figure(figsize=(5, 4))
num_lines = len(experiments)
sns.set_palette(sns.color_palette("coolwarm",num_lines))

ax = plt.axes()
def make_label(experiment):
    return "{:+.1f}, {:+.1f}".format(*experiment)

for experiment in experiments:
    #print experiment
    data = get_data(experiment)
    if data is None:
        print("No data for {}".format(experiment))
        continue
    times = np.array(data.index)
    normalized = data[['CO2']].values / highest_co2
    ax.plot(times, normalized, label=make_label(experiment))
    #normalized.plot.line(ax=ax, label=make_label(experiment))
    #ax.text(times[-10], normalized[-10], make_label(experiment))

#plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.xlabel('Time (s)')
plt.ylabel('Normalized CO$_2$ concentration')
Out[26]:
<matplotlib.text.Text at 0x11ba704d0>
In [27]:
sns.set_palette('Set1')

def my_function(time, rate):
    "Thing we want to fit."
    return 1. - np.exp(-1*time*rate)
import scipy.optimize


def fit_rate(data):
    normalized = data[['CO2']] / highest_co2
    x_data = np.array(normalized.index)
    y_data = normalized.values[:,0] 
    
    try:
        popt, pcov = scipy.optimize.curve_fit(my_function, x_data, y_data)
    except RuntimeError as e:
        print("☝️ couldn't fit: {}".format(e.message))
        return np.nan
        
    optimal_parameters = popt
    parameter_errors = np.sqrt(np.diag(pcov))
    print("Rate: {} +/- {} (1 st. dev.)".format(optimal_parameters[0],parameter_errors[0]))

    plt.plot(x_data, y_data, 'o')
    plt.plot(x_data, my_function(x_data, *optimal_parameters))
    plt.xlabel('Time (s)')
    plt.ylabel('Normalized CO$_2$ concentration')
    
    plt.show()
    fitted_rate = optimal_parameters[0]
    return fitted_rate

rates = []
for experiment in experiments:
    print experiment
    data = get_data(experiment)
    if data is None:
        print("☝️No data")
        rates.append(np.nan)
        continue
    rate = fit_rate(data)
    rates.append(rate)

rates
    
[-7.5 -6.5]
Rate: 0.000504582173939 +/- 3.35252928628e-07 (1 st. dev.)
[-7.5   -5.875]
Rate: 0.150500321688 +/- 0.000945486033578 (1 st. dev.)
[-7.5  -5.25]
Rate: 0.824716324277 +/- 0.0266081078102 (1 st. dev.)
[-7.5   -4.625]
Rate: 1670.63717576 +/- 10.7388828465 (1 st. dev.)
[-7.5 -4. ]
Rate: 4109.52414866 +/- 7.81666120287 (1 st. dev.)
[-7.5   -3.375]
Rate: 4395.8708957 +/- 36.9253473616 (1 st. dev.)
[-7.5  -2.75]
Rate: 35.1917586632 +/- 0.513583486748 (1 st. dev.)
[-7.5   -2.125]
Rate: 1.29560022604 +/- 0.0135245349679 (1 st. dev.)
[-7.5 -1.5]
Rate: 452.479866675 +/- 5.09359117378 (1 st. dev.)
[-6.8125 -6.5   ]
Rate: 0.00347791877856 +/- 1.23293927866e-05 (1 st. dev.)
[-6.8125 -5.875 ]
Rate: 0.131361558227 +/- 0.000877537990693 (1 st. dev.)
[-6.8125 -5.25  ]
Rate: 6.30964506472 +/- 0.111490971072 (1 st. dev.)
[-6.8125 -4.625 ]
Rate: 472.194486022 +/- 2.096863466 (1 st. dev.)
[-6.8125 -4.    ]
Rate: 3122.89211355 +/- 7.52538160841 (1 st. dev.)
[-6.8125 -3.375 ]
Rate: 3969.62713823 +/- 44.8005441464 (1 st. dev.)
[-6.8125 -2.75  ]
Rate: 4.1471420973 +/- 0.0634633990185 (1 st. dev.)
[-6.8125 -2.125 ]
Rate: 0.726029176855 +/- 0.0105147503323 (1 st. dev.)
[-6.8125 -1.5   ]
Rate: 247.130060955 +/- 2.44441218364 (1 st. dev.)
[-6.125 -6.5  ]
Rate: 0.00340388142007 +/- 1.25662625398e-05 (1 st. dev.)
[-6.125 -5.875]
Rate: 0.124968665658 +/- 0.000767722410191 (1 st. dev.)
[-6.125 -5.25 ]
Rate: 1.39300211971 +/- 0.0216902211823 (1 st. dev.)
[-6.125 -4.625]
Rate: 98.9654755498 +/- 0.473364159407 (1 st. dev.)
[-6.125 -4.   ]
Rate: 657.126787343 +/- 3.38823673153 (1 st. dev.)
[-6.125 -3.375]
Rate: 188.870237652 +/- 1.75566554178 (1 st. dev.)
[-6.125 -2.75 ]
Rate: 0.755176859558 +/- 0.0116283186801 (1 st. dev.)
[-6.125 -2.125]
Rate: 0.41837744911 +/- 0.00738437106652 (1 st. dev.)
[-6.125 -1.5  ]
Rate: 262.348850294 +/- 3.01018533348 (1 st. dev.)
[-5.4375 -6.5   ]
Rate: 0.000768155579139 +/- 2.81103633649e-06 (1 st. dev.)
[-5.4375 -5.875 ]
Rate: 0.0094333174238 +/- 3.05541088081e-05 (1 st. dev.)
[-5.4375 -5.25  ]
Rate: 0.0249882643383 +/- 6.97178357252e-05 (1 st. dev.)
[-5.4375 -4.625 ]
Rate: 3.74436109617 +/- 0.0218517122301 (1 st. dev.)
[-5.4375 -4.    ]
Rate: 4.92074292792 +/- 0.0574272447458 (1 st. dev.)
[-5.4375 -3.375 ]
Rate: 93.9868411842 +/- 0.85317175301 (1 st. dev.)
[-5.4375 -2.75  ]
Rate: 0.349837857218 +/- 0.00593033461134 (1 st. dev.)
[-5.4375 -2.125 ]
Rate: 0.197162028292 +/- 0.00283866720558 (1 st. dev.)
[-5.4375 -1.5   ]
Rate: 2.08516470308 +/- 0.0460140029007 (1 st. dev.)
[-4.75 -6.5 ]
Rate: 3.46699483286e-05 +/- 2.96720004417e-07 (1 st. dev.)
[-4.75  -5.875]
Rate: 0.000182275358275 +/- 8.63476456122e-07 (1 st. dev.)
[-4.75 -5.25]
Rate: 0.00015322360536 +/- 7.51722458661e-07 (1 st. dev.)
[-4.75  -4.625]
Rate: 0.114567764435 +/- 0.00108813353474 (1 st. dev.)
[-4.75 -4.  ]
Rate: 0.566122797865 +/- 0.00325521499648 (1 st. dev.)
[-4.75  -3.375]
Rate: 25.4464716246 +/- 0.174325661007 (1 st. dev.)
[-4.75 -2.75]
Rate: 0.0586206258399 +/- 0.000907696849276 (1 st. dev.)
[-4.75  -2.125]
Rate: 0.34574625748 +/- 0.00117678721383 (1 st. dev.)
[-4.75 -1.5 ]
Rate: 1.88817877479 +/- 0.0153086008759 (1 st. dev.)
[-4.0625 -6.5   ]
Rate: 6.12645487706e-07 +/- 5.74499670187e-09 (1 st. dev.)
[-4.0625 -5.875 ]
Rate: 7.72573952358e-07 +/- 5.54280230633e-09 (1 st. dev.)
[-4.0625 -5.25  ]
Rate: 8.85529690587e-07 +/- 7.17188001077e-09 (1 st. dev.)
[-4.0625 -4.625 ]
Rate: 2.37424301963e-06 +/- 1.06463224648e-08 (1 st. dev.)
[-4.0625 -4.    ]
Rate: 0.142550140732 +/- 0.000217562121022 (1 st. dev.)
[-4.0625 -3.375 ]
Rate: 0.0942633906432 +/- 0.000732435228871 (1 st. dev.)
[-4.0625 -2.75  ]
Rate: 0.00303356370505 +/- 1.51110359725e-05 (1 st. dev.)
[-4.0625 -2.125 ]
Rate: 0.0132286148624 +/- 2.13207505955e-05 (1 st. dev.)
[-4.0625 -1.5   ]
Rate: 0.0127463123154 +/- 3.32993558638e-05 (1 st. dev.)
[-3.375 -6.5  ]
☝️ couldn't fit: Optimal parameters not found: Number of calls to function has reached maxfev = 400.
[-3.375 -5.875]
/Users/rwest/anaconda/envs/rmg6/lib/python2.7/site-packages/ipykernel/__main__.py:17: DeprecationWarning: BaseException.message has been deprecated as of Python 2.6
☝️ couldn't fit: Optimal parameters not found: Number of calls to function has reached maxfev = 400.
[-3.375 -5.25 ]
/Users/rwest/anaconda/envs/rmg6/lib/python2.7/site-packages/scipy/optimize/minpack.py:779: OptimizeWarning: Covariance of the parameters could not be estimated
  category=OptimizeWarning)
Rate: -1.1933417947e-10 +/- inf (1 st. dev.)
[-3.375 -4.625]
Rate: -6.80082567044e-09 +/- inf (1 st. dev.)
[-3.375 -4.   ]
Rate: 0.00251870819693 +/- 6.17577520839e-06 (1 st. dev.)
[-3.375 -3.375]
☝️ couldn't fit: Optimal parameters not found: Number of calls to function has reached maxfev = 400.
[-3.375 -2.75 ]
Rate: -6.09315030518e-10 +/- inf (1 st. dev.)
[-3.375 -2.125]
Rate: -1.89653659959e-10 +/- inf (1 st. dev.)
[-3.375 -1.5  ]
☝️ couldn't fit: Optimal parameters not found: Number of calls to function has reached maxfev = 400.
[-2.6875 -6.5   ]
☝️ couldn't fit: Optimal parameters not found: Number of calls to function has reached maxfev = 400.
[-2.6875 -5.875 ]
Rate: -1.26973265672e-10 +/- inf (1 st. dev.)
[-2.6875 -5.25  ]
Rate: -1.41295963179e-10 +/- inf (1 st. dev.)
[-2.6875 -4.625 ]
☝️ couldn't fit: Optimal parameters not found: Number of calls to function has reached maxfev = 400.
[-2.6875 -4.    ]
☝️ couldn't fit: Optimal parameters not found: Number of calls to function has reached maxfev = 400.
[-2.6875 -3.375 ]
☝️ couldn't fit: Optimal parameters not found: Number of calls to function has reached maxfev = 400.
[-2.6875 -2.75  ]
Rate: -3.59205569537e-10 +/- inf (1 st. dev.)
[-2.6875 -2.125 ]
Rate: -8.33977925467e-08 +/- inf (1 st. dev.)
[-2.6875 -1.5   ]
☝️ couldn't fit: Optimal parameters not found: Number of calls to function has reached maxfev = 400.
[-2.  -6.5]
Rate: -2.85116384348e-10 +/- inf (1 st. dev.)
[-2.    -5.875]
☝️ couldn't fit: Optimal parameters not found: Number of calls to function has reached maxfev = 400.
[-2.   -5.25]
Rate: -1.26076967227e-10 +/- inf (1 st. dev.)
[-2.    -4.625]
☝️ couldn't fit: Optimal parameters not found: Number of calls to function has reached maxfev = 400.
[-2. -4.]
Rate: -1.96144793303e-10 +/- inf (1 st. dev.)
[-2.    -3.375]
☝️ couldn't fit: Optimal parameters not found: Number of calls to function has reached maxfev = 400.
[-2.   -2.75]
Rate: -1.00113114352e-09 +/- inf (1 st. dev.)
[-2.    -2.125]
Rate: -2.17361630741e-07 +/- inf (1 st. dev.)
[-2.  -1.5]
Rate: -8.00675265678e-10 +/- inf (1 st. dev.)
Out[27]:
[0.00050458217393925162,
 0.15050032168801367,
 0.82471632427714237,
 1670.6371757618408,
 4109.5241486580717,
 4395.8708956951577,
 35.191758663181695,
 1.295600226043585,
 452.4798666749324,
 0.0034779187785619137,
 0.13136155822666262,
 6.3096450647213871,
 472.19448602242409,
 3122.8921135505452,
 3969.6271382267341,
 4.1471420972993238,
 0.72602917685490442,
 247.13006095491411,
 0.0034038814200712583,
 0.12496866565792325,
 1.393002119712601,
 98.965475549761095,
 657.12678734335327,
 188.87023765184887,
 0.75517685955811742,
 0.4183774491101302,
 262.34885029401158,
 0.00076815557913949241,
 0.0094333174238047015,
 0.024988264338293154,
 3.7443610961748548,
 4.9207429279169501,
 93.986841184219287,
 0.34983785721830341,
 0.19716202829242518,
 2.0851647030845837,
 3.4669948328584798e-05,
 0.00018227535827506226,
 0.00015322360536028075,
 0.11456776443525636,
 0.56612279786499675,
 25.446471624550991,
 0.058620625839867088,
 0.34574625748009169,
 1.888178774794125,
 6.1264548770627916e-07,
 7.725739523577382e-07,
 8.8552969058738659e-07,
 2.3742430196335243e-06,
 0.14255014073150796,
 0.094263390643230027,
 0.0030335637050512682,
 0.013228614862435395,
 0.012746312315404016,
 nan,
 nan,
 -1.1933417947007279e-10,
 -6.8008256704396179e-09,
 0.0025187081969282972,
 nan,
 -6.0931503051839586e-10,
 -1.8965365995908334e-10,
 nan,
 nan,
 -1.2697326567212617e-10,
 -1.4129596317867359e-10,
 nan,
 nan,
 nan,
 -3.5920556953721034e-10,
 -8.3397792546709147e-08,
 nan,
 -2.8511638434809939e-10,
 nan,
 -1.2607696722716903e-10,
 nan,
 -1.961447933033468e-10,
 nan,
 -1.0011311435190012e-09,
 -2.1736163074128254e-07,
 -8.0067526567831072e-10]
In [28]:
rates = np.array(rates)
fixed_rates = rates * (rates>0) + (1e-9 * (rates<0))
log_rates = np.log(fixed_rates)
log_rates
/Users/rwest/anaconda/envs/rmg6/lib/python2.7/site-packages/ipykernel/__main__.py:2: RuntimeWarning: invalid value encountered in greater
  from ipykernel import kernelapp as app
/Users/rwest/anaconda/envs/rmg6/lib/python2.7/site-packages/ipykernel/__main__.py:2: RuntimeWarning: invalid value encountered in less
  from ipykernel import kernelapp as app
Out[28]:
array([ -7.59177985,  -1.89379006,  -0.1927158 ,   7.42096038,
         8.32106252,   8.38842095,   3.56081193,   0.25897408,
         6.11474327,  -5.66132122,  -2.02980177,   1.84207943,
         6.15739095,   8.04651481,   8.28642745,   1.42241945,
        -0.32016508,   5.50991476,  -5.6828389 ,  -2.07969225,
         0.33146122,   4.59477106,   6.48787698,   5.24106021,
        -0.28080331,  -0.87137127,   5.56967511,  -7.17151827,
        -4.66350745,  -3.68934899,   1.320251  ,   1.59345952,
         4.54315479,  -1.0502855 ,  -1.62372941,   0.73484785,
       -10.26963729,  -8.60999206,  -8.78361223,  -2.1665888 ,
        -0.56894427,   3.23657709,  -2.83666867,  -1.06205013,
         0.63561275, -14.30547939, -14.0735381 , -13.93707985,
       -12.9508319 ,  -1.94806148,  -2.36166239,  -5.79801721,
        -4.325373  ,  -4.36251328,          nan,          nan,
       -20.72326584, -20.72326584,  -5.98400913,          nan,
       -20.72326584, -20.72326584,          nan,          nan,
       -20.72326584, -20.72326584,          nan,          nan,
                nan, -20.72326584, -20.72326584,          nan,
       -20.72326584,          nan, -20.72326584,          nan,
       -20.72326584,          nan, -20.72326584, -20.72326584, -20.72326584])
In [29]:
rate_grid = np.reshape(log_rates, (grid_size,grid_size))
In [30]:
ax = sns.heatmap(rate_grid)
In [31]:
extent = carbon_range + oxygen_range
extent
Out[31]:
(-7.5, -2, -6.5, -1.5)
In [32]:
# Because the center of a corner pixel is in fact the corner of the grid
# we want to stretch the image a little
c_step = mesh[0,1,0]-mesh[0,0,0]
o_step = mesh[1,0,1]-mesh[1,0,0]
carbon_range2 = (carbon_range[0]-c_step/2, carbon_range[1]+c_step/2)
oxygen_range2 = (oxygen_range[0]-c_step/2, oxygen_range[1]+c_step/2)
extent2 = carbon_range2 + oxygen_range2
extent2
Out[32]:
(-7.84375, -1.65625, -6.84375, -1.15625)
In [33]:
plt.imshow(rate_grid.T, interpolation='none', origin='lower', extent=extent2, aspect='equal')
plt.plot(-5.997, -4.485, 'ok')
plt.text(-5.997, -4.485, 'Ni(111)')
Out[33]:
<matplotlib.text.Text at 0x11c32af50>
In [34]:
# Binding energies extracted from:
u"""
Medford, A. J.; Lausche, A. C.; Abild-Pedersen, F.;
Temel, B.; Schjødt, N. C.; Nørskov, J. K.; Studt, F. 
Activity and Selectivity Trends in 
Synthesis Gas Conversion to Higher Alcohols. 
Topics in Catalysis 2014, 57 (1-4), 135–142 
DOI: 10.1007/s11244-013-0169-0.
"""
medford_energies = { # Carbon, then Oxygen
'Ru': ( 0.010349288486416697, -2.8153856448231256),
'Rh': ( 0.16558861578266493, -2.546620868091181),
'Ni': ( 0.3001293661060802, -2.5881741535441853),
'Ir': ( 0.36222509702457995, -2.826185484230718),
'Pd': ( 0.28460543337645516, -1.207119596734621),
'Pt': ( 0.8796895213454077, -1.445820136503547),
'Cu': ( 2.323415265200518, -1.7218249542757729),
'Ag': ( 3.855109961190168, -0.8341504215550701),
'Au': ( 3.5601552393272975, -0.10963108355266138),
}
In [35]:
# Shift Medford's energies so that Ni matches Wayne Blaylock's Ni
blaylock_ni = np.array([-5.997, -4.485])
old_ni = np.array(medford_energies['Ni'])
shifted_energies = {metal: tuple(blaylock_ni + np.array(E)-old_ni) for metal,E in medford_energies.items()}
shifted_energies
Out[35]:
{'Ag': (-2.4420194049159121, -2.7309762680108851),
 'Au': (-2.7369741267787826, -2.006456930008476),
 'Cu': (-3.9737141009055619, -3.6186508007315878),
 'Ir': (-5.9349042690814997, -4.723011330686532),
 'Ni': (-5.9969999999999999, -4.4849999999999994),
 'Pd': (-6.0125239327296249, -3.1039454431904363),
 'Pt': (-5.417439844760672, -3.3426459829593616),
 'Rh': (-6.1315407503234152, -4.4434467145469956),
 'Ru': (-6.2867800776196638, -4.712211491278941)}
In [36]:
plt.imshow(rate_grid.T, interpolation='none', origin='lower', extent=extent2, aspect='equal')
for metal, coords in shifted_energies.iteritems():
    plt.plot(coords[0], coords[1], 'ok')
    plt.text(coords[0], coords[1], metal)
plt.xlim(carbon_range)
plt.ylim(oxygen_range)
Out[36]:
(-6.5, -1.5)
In [37]:
# Binding energies for close packed surfaces from
"""
Abild-Pedersen, F.; Greeley, J.; Studt, F.; Rossmeisl, J.;
Munter, T. R.; Moses, P. G.; Skúlason, E.; Bligaard, T.; Norskov, J. K. 
Scaling Properties of Adsorption Energies for Hydrogen-Containing 
Molecules on Transition-Metal Surfaces. 
Phys. Rev. Lett. 2007, 99 (1), 016105 
DOI: 10.1103/PhysRevLett.99.016105.
"""
abildpedersen_energies = { # Carbon, then Oxygen
'Ru': ( -6.397727272727272, -5.104763568600047),
'Rh': ( -6.5681818181818175, -4.609771721406942),
'Ni': ( -6.045454545454545, -4.711681807593758),
'Ir': ( -6.613636363636363, -5.94916142557652),
'Pd': ( -6, -3.517877940833916),
'Pt': ( -6.363636363636363, -3.481481481481482),
'Cu': ( -4.159090909090907, -3.85272536687631),
'Ag': ( -2.9545454545454533, -2.9282552993244817),
'Au': ( -3.7499999999999973, -2.302236198462614),
}
In [38]:
abildpedersen_energies['Pt'][0] - abildpedersen_energies['Ni'][0]
Out[38]:
-0.31818181818181834
In [39]:
plt.imshow(rate_grid.T, interpolation='none', origin='lower', extent=extent2, aspect='equal')
for metal, coords in abildpedersen_energies.iteritems():
    color = {'Ag':'w','Au':'w','Cu':'w'}.get(metal,'k')
    plt.plot(coords[0], coords[1], 'o'+color)
    plt.text(coords[0], coords[1], metal, color=color)
plt.xlim(carbon_range)
plt.ylim(oxygen_range)
plt.xlabel('$\Delta E^C$ (eV)')
plt.ylabel('$\Delta E^O$ (eV)')
Out[39]:
<matplotlib.text.Text at 0x11d730150>

Alternative rate definition

For this, we plot the CO2 concentration profiles on a log-log plot, and see the characteristic time as the time at which it crosses the diagonal, i.e. when does $\left( \log_{10}(time) + \log_{10}(\frac{CO_2}{CO_2max}) \right) \ge -10$

In [40]:
"""
Plot everything on a log-log plot
"""
import seaborn as sns
plt.figure(figsize=(5, 4))
num_lines = len(experiments)
sns.set_palette(sns.color_palette("coolwarm",num_lines))
ax = plt.axes()

def make_label(experiment):
    return "{:+.1f}, {:+.1f}".format(*experiment)

for experiment in experiments:
    print experiment
    data = get_data(experiment)
    if data is None:
        print("☝️No data")
        continue
    times = np.array(data.index)
    normalized = data[['CO2']].values[:,0] / highest_co2
    ax.loglog(times, normalized, label=make_label(experiment))
    try:
        i = (np.nonzero((np.log10(times)+np.log10(normalized)) > -10))[0][0]
    except IndexError:
        i = -1
    print i, times[i], normalized[i]
    ax.text(times[i], normalized[i], make_label(experiment), rotation=45, ha='center', va='center', fontsize=12)
    
plt.plot([1e-10, 1],[1,1e-10], 'g:')
plt.ylim(1e-10, 1)
plt.xlim(1e-10, 1)
#plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.xlabel('Time (s)')
plt.ylabel('Normalized CO2 concentration')
[-7.5 -6.5]
5090 0.000513578082952 1.94819347457e-07
[-7.5   -5.875]
/Users/rwest/anaconda/envs/rmg6/lib/python2.7/site-packages/ipykernel/__main__.py:23: RuntimeWarning: divide by zero encountered in log10
4565 1.85863315708e-05 5.39774385091e-06
[-7.5  -5.25]
3093 1.3209964688e-06 7.61098920181e-05
[-7.5   -4.625]
2855 1.06937134557e-06 9.36913399659e-05
[-7.5 -4. ]
2884 7.91201429063e-07 0.000127213165266
[-7.5   -3.375]
2622 2.35510510273e-07 0.000430380911287
[-7.5  -2.75]
3113 1.85919034089e-07 0.000538613068765
[-7.5   -2.125]
2717 1.98501169593e-07 0.000507799589196
[-7.5 -1.5]
2879 4.11672091604e-07 0.000242936100676
[-6.8125 -6.5   ]
4752 0.000122577717674 8.20796242217e-07
[-6.8125 -5.875 ]
2911 5.26549829718e-06 1.91126086508e-05
[-6.8125 -5.25  ]
2916 1.58962283825e-06 6.33820632805e-05
[-6.8125 -4.625 ]
3076 1.80373752033e-06 5.58383676836e-05
[-6.8125 -4.    ]
2773 9.35478084833e-07 0.000107252734655
[-6.8125 -3.375 ]
3012 2.97842489699e-07 0.000339916267995
[-6.8125 -2.75  ]
2612 2.76273549078e-07 0.00036665289471
[-6.8125 -2.125 ]
2951 2.98509448669e-07 0.000335496808129
[-6.8125 -1.5   ]
2505 1.1923824773e-06 8.38756938659e-05
[-6.125 -6.5  ]
4315 7.55485504653e-05 1.32909943718e-06
[-6.125 -5.875]
2883 5.62722797126e-06 1.77879382187e-05
[-6.125 -5.25 ]
3246 4.2859265821e-06 2.34282205967e-05
[-6.125 -4.625]
3466 4.72077539751e-06 2.12535566181e-05
[-6.125 -4.   ]
3332 1.46582160503e-06 6.87047592332e-05
[-6.125 -3.375]
2866 4.9034920186e-07 0.000204279105215
[-6.125 -2.75 ]
2988 4.70187706384e-07 0.000216032579675
[-6.125 -2.125]
3021 6.32103811824e-07 0.000158348258477
[-6.125 -1.5  ]
3008 4.82611039926e-06 2.09875832195e-05
[-5.4375 -6.5   ]
4760 9.39072550941e-05 1.06669056498e-06
[-5.4375 -5.875 ]
3348 1.41805353639e-05 7.07345503921e-06
[-5.4375 -5.25  ]
3306 1.58923646745e-05 6.31056807501e-06
[-5.4375 -4.625 ]
3928 1.29752385841e-05 7.72975133471e-06
[-5.4375 -4.    ]
3603 2.50233654835e-06 4.00997901998e-05
[-5.4375 -3.375 ]
3158 8.251863997e-07 0.000122372193089
[-5.4375 -2.75  ]
2791 1.01801528936e-06 9.85068524764e-05
[-5.4375 -2.125 ]
2708 2.58963637826e-06 3.88200334805e-05
[-5.4375 -1.5   ]
3963 1.93120382267e-05 5.2188518303e-06
[-4.75 -6.5 ]
4228 0.000141735429518 7.0609679714e-07
[-4.75  -5.875]
3796 6.83380212181e-05 1.46373511855e-06
[-4.75 -5.25]
4426 8.90920637461e-05 1.13007895187e-06
[-4.75  -4.625]
4658 3.9367808746e-05 2.55028830698e-06
[-4.75 -4.  ]
3702 5.9635128114e-06 1.68428118558e-05
[-4.75  -3.375]
3959 2.07704526155e-06 4.84884935656e-05
[-4.75 -2.75]
3988 3.80934067453e-06 2.63806162762e-05
[-4.75  -2.125]
3190 3.529659005e-05 2.83413251452e-06
[-4.75 -1.5 ]
3561 8.90183895077e-05 1.13923935584e-06
[-4.0625 -6.5   ]
5009 0.000992626484175 1.01098644663e-07
[-4.0625 -5.875 ]
4582 0.000944571523009 1.06373325801e-07
[-4.0625 -5.25  ]
5630 0.00118682619781 8.4327431088e-08
[-4.0625 -4.625 ]
5599 0.00112129787063 8.958333406e-08
[-4.0625 -4.    ]
4181 3.63925014801e-05 2.75324425396e-06
[-4.0625 -3.375 ]
3809 1.06810475933e-05 9.36296621109e-06
[-4.0625 -2.75  ]
3775 2.60822611352e-05 3.83903118916e-06
[-4.0625 -2.125 ]
4407 0.000103343169595 9.68699467621e-07
[-4.0625 -1.5   ]
3818 0.00116810978596 8.66579683672e-08
[-3.375 -6.5  ]
-1 1.00193260469 0.0
[-3.375 -5.875]
-1 1.00090445215 0.0
[-3.375 -5.25 ]
-1 1.00028638516 2.43520568266e-11
[-3.375 -4.625]
-1 0.018720391576 2.71928531984e-10
[-3.375 -4.   ]
5264 0.00076669061667 1.31065534352e-07
[-3.375 -3.375]
-1 0.00144386513779 0.0
[-3.375 -2.75 ]
-1 1.01245933502 0.0
[-3.375 -2.125]
-1 1.00055160153 8.51297031636e-12
[-3.375 -1.5  ]
-1 0.0354807170287 4.20187775724e-16
[-2.6875 -6.5   ]
-1 1.00066715114 0.0
[-2.6875 -5.875 ]
-1 1.00216684184 0.0
[-2.6875 -5.25  ]
-1 1.00118297556 0.0
[-2.6875 -4.625 ]
-1 1.00024160783 0.0
[-2.6875 -4.    ]
-1 1.00100257514 0.0
[-2.6875 -3.375 ]
-1 1.00159553428 0.0
[-2.6875 -2.75  ]
-1 1.00479395746 0.0
[-2.6875 -2.125 ]
-1 0.001845890325 0.0
[-2.6875 -1.5   ]
-1 0.263938204963 0.0
[-2.  -6.5]
-1 1.00096093055 0.0
[-2.    -5.875]
-1 1.00065531449 0.0
[-2.   -5.25]
-1 1.00072001947 0.0
[-2.    -4.625]
-1 1.00087141269 0.0
[-2. -4.]
-1 1.00018999866 0.0
[-2.    -3.375]
-1 1.00016546685 0.0
[-2.   -2.75]
-1 1.00351642812 0.0
[-2.    -2.125]
-1 0.000715136208279 0.0
[-2.  -1.5]
-1 0.179084350044 0.0
Out[40]:
<matplotlib.text.Text at 0x120f73b50>
In [53]:
"""
An alternative way of defining rates.
"""
new_rates = []
for experiment in experiments:
    #print experiment
    data = get_data(experiment)
    if data is None:
        print("No data for experiment {}".format(experiment))
        new_rates.append(np.nan)
        continue
        
    times = np.array(data.index)
    normalized = data[['CO2']].values[:,0] / highest_co2
    try:
        i = (np.nonzero((np.log10(times)+np.log10(normalized)) > -10))[0][0]
        time = times[i]
    except IndexError:
        time = 1.0
    new_rates.append(1./time)
new_log_rates = np.log(np.array(new_rates))
print new_log_rates
/Users/rwest/anaconda/envs/rmg6/lib/python2.7/site-packages/ipykernel/__main__.py:16: RuntimeWarning: divide by zero encountered in log10
[  8.88812052  12.21029427  14.76270805  14.87947389  15.08624509
  16.21831815  16.84921209  16.64456016  15.81949745  10.71542947
  13.65752799  14.48731689  14.3883084   14.90379989  16.06640208
  16.57343441  16.43058327  14.69063118  11.61602318  13.55333138
  13.68865812  13.62600433  14.50521814  15.5483624   15.93893469
  15.57257806  13.17091627  11.53318412  12.82398877  12.71029781
  12.84738723  14.12747637  15.05902059  15.07335293  13.94687101
  11.82996751  11.19205501  11.65517124  11.46507434  12.06800942
  13.60419676  14.4122212   13.50562998  11.18038971  10.43106783
   9.75224944   9.78089457   9.59614895   9.5204914   12.47324586
  13.17486708  11.69462433  10.29081576   7.57358478   0.           0.
   2.66247466   5.06728078  10.13993801   0.           0.           0.79207734
   0.           0.           0.           0.           0.           0.           0.
   0.           0.           0.           0.           0.           0.           0.
   0.           0.           0.           0.           0.        ]
In [54]:
new_rate_grid = np.reshape(new_log_rates, (grid_size,grid_size))
plt.imshow(new_rate_grid.T, interpolation='none', origin='lower', extent=extent2, aspect='equal')
for metal, coords in abildpedersen_energies.iteritems():
    color = {'Ag':'w','Au':'w',}.get(metal,'k')
    plt.plot(coords[0], coords[1], 'o'+color)
    plt.text(coords[0], coords[1], metal, color=color)
plt.xlim(carbon_range)
plt.ylim(oxygen_range)
plt.xlabel('$\Delta E^C$ (eV)')
plt.ylabel('$\Delta E^O$ (eV)')
Out[54]:
<matplotlib.text.Text at 0x119af2dd0>
In [55]:
new_rate_grid = np.reshape(new_log_rates, (grid_size,grid_size))
plt.imshow(new_rate_grid.T, interpolation='spline16', origin='lower', extent=extent2, aspect='equal')
for metal, coords in abildpedersen_energies.iteritems():
    color = {'Ag':'w','Au':'w',}.get(metal,'k')
    plt.plot(coords[0], coords[1], 'o'+color)
    plt.text(coords[0], coords[1], metal, color=color)
plt.xlim(carbon_range)
plt.ylim(oxygen_range)
plt.xlabel('$\Delta E^C$ (eV)')
plt.ylabel('$\Delta E^O$ (eV)')
Out[55]:
<matplotlib.text.Text at 0x14e0239d0>

Count the reactions

In [43]:
def get_reaction_count(experiment):
    d = directory(*experiment)
    f = os.path.join(d,'chemkin','chem_annotated.inp')
    with open(f) as chemkin:
        r = re.compile('\! Reaction index: Chemkin \#(\d+)')
        for line in chemkin:
            m = r.match(line)
            if m:
                count = int(m.group(1))
    return count
    
get_reaction_count(experiment)
Out[43]:
4059
In [44]:
i=35
print experiments[i]
print get_reaction_count(experiments[i])
[-5.4375 -1.5   ]
534
In [45]:
reaction_counts = map(get_reaction_count, experiments)
reaction_counts_grid = np.log10(np.reshape(reaction_counts, (grid_size,grid_size)))

plt.imshow(reaction_counts_grid.T, interpolation='none', origin='lower', extent=extent2, aspect='equal')
for metal, coords in abildpedersen_energies.iteritems():
    color = {'Ag':'w','Au':'w','Cu':'w'}.get(metal,'k')
    color='k'
    plt.plot(coords[0], coords[1], 'o'+color)
    plt.text(coords[0], coords[1], metal, color=color)
plt.xlim(carbon_range2)
plt.ylim(oxygen_range2)
plt.xlabel('$\Delta E^C$ (eV)')
plt.ylabel('$\Delta E^O$ (eV)')

for e,n in zip(experiments,reaction_counts):
    plt.text(e[0],e[1],n,color='w',ha='center', va='center')
#plt.colorbar()
In [46]:
# A linear one, just to check it looks the same
reaction_counts_grid = np.reshape(reaction_counts, (grid_size,grid_size))
ax = sns.heatmap(reaction_counts_grid.T[::-1,:], annot=True, fmt='d', square=True)